#include "glni.h"


double glni::interpolate1d_linear(double x1, double x2, double y1, double y2, double x)
{
	if (x1 == x2) return 0.5*(y1 + y2);
	else return (x - x1)*(y2 -y1)/(x2 - x1) + y1;
}

glni::line::line()
{
	initialized_ = false;
	size_ = -1;
	nodes_ = nullptr;
	weights_ = nullptr;
}

glni::line::line(unsigned int size) : line()
{
	initiate(size);
}

glni::line::~line()
{
	clear();
}

void glni::line::clear()
{
	if (nodes_ != nullptr) delete[] nodes_;
	if (weights_ != nullptr) delete[] weights_;
	initialized_ = false;
	size_ = -1;
	nodes_ = nullptr;
	weights_ = nullptr;
	return;
}

double glni::line::node(unsigned int idx)
{
#ifndef GLNI_NO_CHECK
	if (!initialized_)
		throw std::invalid_argument("glni::line is not initialized.");

	if (idx >= size_)
		throw std::invalid_argument("glni::line::node(...) got an invalid index number.");
#endif

	return nodes_[idx];
}

double glni::line::weight(unsigned int idx)
{
#ifndef GLNI_NO_CHECK
	if (!initialized_)
		throw std::invalid_argument("glni::line is not initialized.");

	if (idx >= size_)
		throw std::invalid_argument("glni::line::weight(...) got an invalid index number.");
#endif

	return weights_[idx];
}

double glni::line::integral(double low, double hig, func_1d_ptr fx, void *att, void *instance)
{
#ifndef GLNI_NO_CHECK
	if (!initialized_)
		throw std::invalid_argument("glni::line is not initialized.");

	if (hig <= low)
		throw std::invalid_argument("glni::line::integral(...) got an invalid integration range.");
#endif

	double scl = (hig - low)/2.0;
	double gx, inte = 0.0;
	for (int i = 0; i < size_; i++)
	{
		gx = 0.5*(low + hig) + scl*nodes_[i];
		inte += weights_[i]*scl*fx(gx, low, hig, att, instance);
	}
	return inte;
}

double glni::line::integral(const std::vector<double> &ix, const std::vector<double> &iv, func_1d_ptr2 fx, void *att, void *instance)
{
#ifndef GLNI_NO_CHECK
	if (!initialized_)
		throw std::invalid_argument("glni::line is not initialized.");

	if (ix.size() <= 1 || iv.size() <= 1 || ix.size() != iv.size())
		throw std::invalid_argument("glni::line got an invalid vector.");

	for (int i = 0; i < ix.size()-1; i++)
	{
		if (ix[i] >= ix[i+1])
			throw std::runtime_error("glni::line::integral(...) got an invalid integration range.");
	}
#endif

	double scl = (ix.back() - ix.front())/2.0;
	double gx, inte = 0.0;
	for (int i = 0; i < size_; i++)
	{
		gx = 0.5*(ix.front() + ix.back()) + scl*nodes_[i];
		inte += weights_[i]*scl*fx(gx, ix, iv, att, instance);
	}
	return inte;
}

void glni::line::initiate(unsigned int size)
{
#ifndef GLNI_NO_CHECK
	if (size == 0 || size > 20)
		throw std::invalid_argument("glni::line::initiate(...) got an invalid initiating order.");
#endif

	initialized_ = true;
	size_ = size;
	nodes_ = new double [size_];
	weights_ = new double [size_];

	switch (size_)
	{
	case 1:
		nodes_[0] = 0.0;
		weights_[0] = 2.0;
		break;
	case 2:
		nodes_[0] = -0.577350269189625764509148780502; nodes_[1] = -1.0*nodes_[0];
		weights_[0] = weights_[1] = 1.0;
		break;
	case 3:
		nodes_[0] = -0.774596669241483377035853079956; nodes_[2] = -1.0*nodes_[0];
		nodes_[1] = 0.0;
		weights_[0] = weights_[2] = 0.5555555555555555555555555555565;
		weights_[1] = 0.8888888888888888888888888888889;
		break;
	case 4:
		nodes_[0] = -0.861136311594052575223946488893; nodes_[3] = -1.0*nodes_[0];
		nodes_[1] = -0.339981043584856264802665759103; nodes_[2] = -1.0*nodes_[1];
		weights_[0] = weights_[3] = 0.347854845137453857373063949222;
		weights_[1] = weights_[2] = 0.652145154862546142626936050778;
		break;
	case 5:
		nodes_[0] = -0.906179845938663992797626878299; nodes_[4] = -1.0*nodes_[0];
		nodes_[1] = -0.538469310105683091036314420700; nodes_[3] = -1.0*nodes_[1];
		nodes_[2] = 0.0;
		weights_[0] = weights_[4] = 0.236926885056189087514264040720;
		weights_[1] = weights_[3] = 0.478628670499366468041291514836;
		weights_[2] = 0.568888888888888888888888888889;
		break;
	case 6:
		nodes_[0] = -0.932469514203152027812301554494; nodes_[5] = -1.0*nodes_[0];
		nodes_[1] = -0.661209386466264513661399595020; nodes_[4] = -1.0*nodes_[1];
		nodes_[2] = -0.238619186083196908630501721681; nodes_[3] = -1.0*nodes_[2];
		weights_[0] = weights_[5] = 0.171324492379170345040296142173;
		weights_[1] = weights_[4] = 0.360761573048138607569833513838;
		weights_[2] = weights_[3] = 0.467913934572691047389870343990;
		break;
	case 7:
		nodes_[0] = -0.949107912342758524526189684048; nodes_[6] = -1.0*nodes_[0];
		nodes_[1] = -0.741531185599394439863864773281; nodes_[5] = -1.0*nodes_[1];
		nodes_[2] = -0.405845151377397166906606412077; nodes_[4] = -1.0*nodes_[2];
		nodes_[3] = 0.0;
		weights_[0] = weights_[6] = 0.129484966168869693270611432679;
		weights_[1] = weights_[5] = 0.279705391489276667901467771424;
		weights_[2] = weights_[4] = 0.381830050505118944950369775489;
		weights_[3] = 0.417959183673469387755102040816;
		break;
	case 8:
		nodes_[0] = -0.960289856497536231683560868569; nodes_[7] = -1.0*nodes_[0];
		nodes_[1] = -0.796666477413626739591553936476; nodes_[6] = -1.0*nodes_[1];
		nodes_[2] = -0.525532409916328985817739049189; nodes_[5] = -1.0*nodes_[2];
		nodes_[3] = -0.183434642495649804939476142360; nodes_[4] = -1.0*nodes_[3];
		weights_[0] = weights_[7] = 0.101228536290376259152531354310;
		weights_[1] = weights_[6] = 0.222381034453374470544355994426;
		weights_[2] = weights_[5] = 0.313706645877887287337962201987;
		weights_[3] = weights_[4] = 0.362683783378361982965150449277;
		break;
	case 9:
		nodes_[0] = -0.968160239507626089835576202904; nodes_[8] = -1.0*nodes_[0];
		nodes_[1] = -0.836031107326635794299429788070; nodes_[7] = -1.0*nodes_[1];
		nodes_[2] = -0.613371432700590397308702039341; nodes_[6] = -1.0*nodes_[2];
		nodes_[3] = -0.324253423403808929038538014643; nodes_[5] = -1.0*nodes_[3];
		nodes_[4] = 0.0;
		weights_[0] = weights_[8] = 0.0812743883615744119718921581105;
		weights_[1] = weights_[7] = 0.180648160694857404058472031243;
		weights_[2] = weights_[6] = 0.260610696402935462318742869419;
		weights_[3] = weights_[5] = 0.312347077040002840068630406584;
		weights_[4] = 0.330239355001259763164525069287;
		break;
	case 10:
		nodes_[0] = -0.973906528517171720077964012084; nodes_[9] = -1.0*nodes_[0];
		nodes_[1] = -0.865063366688984510732096688423; nodes_[8] = -1.0*nodes_[1];
		nodes_[2] = -0.679409568299024406234327365115; nodes_[7] = -1.0*nodes_[2];
		nodes_[3] = -0.433395394129247290799265943166; nodes_[6] = -1.0*nodes_[3];
		nodes_[4] = -0.148874338981631210884826001130; nodes_[5] = -1.0*nodes_[4];
		weights_[0] = weights_[9] = 0.0666713443086881375935688098933;
		weights_[1] = weights_[8] = 0.149451349150580593145776339658;
		weights_[2] = weights_[7] = 0.219086362515982043995534934228;
		weights_[3] = weights_[6] = 0.269266719309996355091226921569;
		weights_[4] = weights_[5] = 0.295524224714752870173892994651;
		break;
	case 11:
		nodes_[0]=-0.978228658146056992803938001123; nodes_[10] = -1.0*nodes_[0];
		nodes_[1]=-0.887062599768095299075157769304; nodes_[9] = -1.0*nodes_[1];
		nodes_[2]=-0.730152005574049324093416252031; nodes_[8] = -1.0*nodes_[2];
		nodes_[3]=-0.519096129206811815925725669459; nodes_[7] = -1.0*nodes_[3];
		nodes_[4]=-0.269543155952344972331531985401; nodes_[6] = -1.0*nodes_[4];
		nodes_[5]= 0.0;
		weights_[0]= weights_[10] = 0.0556685671161736664827537204425;
		weights_[1]= weights_[9] = 0.125580369464904624634694299224;
		weights_[2]= weights_[8] = 0.186290210927734251426097641432;
		weights_[3]= weights_[7] = 0.233193764591990479918523704843;
		weights_[4]= weights_[6] = 0.262804544510246662180688869891;
		weights_[5]= 0.272925086777900630714483528336;
		break;
	case 12:
		nodes_[0]=-0.981560634246719250690549090149; nodes_[11] = -1.0*nodes_[0];
		nodes_[1]=-0.904117256370474856678465866119; nodes_[10] = -1.0*nodes_[1];
		nodes_[2]=-0.769902674194304687036893833213; nodes_[9] = -1.0*nodes_[2];
		nodes_[3]=-0.587317954286617447296702418941; nodes_[8] = -1.0*nodes_[3];
		nodes_[4]=-0.367831498998180193752691536644; nodes_[7] = -1.0*nodes_[4];
		nodes_[5]=-0.125233408511468915472441369464; nodes_[6] = -1.0*nodes_[5];
		weights_[0]= weights_[11] = 0.0471753363865118271946159614850;
		weights_[1]= weights_[10] = 0.106939325995318430960254718194;
		weights_[2]= weights_[9] = 0.160078328543346226334652529543;
		weights_[3]= weights_[8] = 0.203167426723065921749064455810;
		weights_[4]= weights_[7] = 0.233492536538354808760849898925;
		weights_[5]= weights_[6] = 0.249147045813402785000562436043;
		break;
	case 13:
		nodes_[0]=-0.984183054718588149472829448807; nodes_[12] = -1.0*nodes_[0];
		nodes_[1]=-0.917598399222977965206547836501; nodes_[11] = -1.0*nodes_[1];
		nodes_[2]=-0.801578090733309912794206489583; nodes_[10] = -1.0*nodes_[2];
		nodes_[3]=-0.642349339440340220643984606996; nodes_[9] = -1.0*nodes_[3];
		nodes_[4]=-0.448492751036446852877912852128; nodes_[8] = -1.0*nodes_[4];
		nodes_[5]=-0.230458315955134794065528121098; nodes_[7] = -1.0*nodes_[5];
		nodes_[6]= 0.0;
		weights_[0]= weights_[12] = 0.0404840047653158795200215922010;
		weights_[1]= weights_[11] = 0.0921214998377284479144217759538;
		weights_[2]= weights_[10] = 0.138873510219787238463601776869;
		weights_[3]= weights_[9] = 0.178145980761945738280046691996;
		weights_[4]= weights_[8] = 0.207816047536888502312523219306;
		weights_[5]= weights_[7] = 0.226283180262897238412090186040;
		weights_[6]= 0.232551553230873910194589515269;
		break;
	case 14:
		nodes_[0]=-0.986283808696812338841597266704; nodes_[13] = -1.0*nodes_[0];
		nodes_[1]=-0.928434883663573517336391139378; nodes_[12] = -1.0*nodes_[1];
		nodes_[2]=-0.827201315069764993189794742650; nodes_[11] = -1.0*nodes_[2];
		nodes_[3]=-0.687292904811685470148019803019; nodes_[10] = -1.0*nodes_[3];
		nodes_[4]=-0.515248636358154091965290718551; nodes_[9] = -1.0*nodes_[4];
		nodes_[5]=-0.319112368927889760435671824168; nodes_[8] = -1.0*nodes_[5];
		nodes_[6]=-0.108054948707343662066244650220; nodes_[7] = -1.0*nodes_[6];
		weights_[0]= weights_[13] = 0.0351194603317518630318328761382;
		weights_[1]= weights_[12] = 0.0801580871597602098056332770629;
		weights_[2]= weights_[11] = 0.121518570687903184689414809072;
		weights_[3]= weights_[10] = 0.157203167158193534569601938624;
		weights_[4]= weights_[9] = 0.185538397477937813741716590125;
		weights_[5]= weights_[8] = 0.205198463721295603965924065661;
		weights_[6]= weights_[7] = 0.215263853463157790195876443316;
		break;
	case 15:
		nodes_[0]=-0.987992518020485428489565718587; nodes_[14] = -1.0*nodes_[0];
		nodes_[1]=-0.937273392400705904307758947710; nodes_[13] = -1.0*nodes_[1];
		nodes_[2]=-0.848206583410427216200648320774; nodes_[12] = -1.0*nodes_[2];
		nodes_[3]=-0.724417731360170047416186054614; nodes_[11] = -1.0*nodes_[3];
		nodes_[4]=-0.570972172608538847537226737254; nodes_[10] = -1.0*nodes_[4];
		nodes_[5]=-0.394151347077563369897207370981; nodes_[9] = -1.0*nodes_[5];
		nodes_[6]=-0.201194093997434522300628303395; nodes_[8] = -1.0*nodes_[6];
		nodes_[7]= 0.0;
		weights_[0]= weights_[14] = 0.0307532419961172683546283935772;
		weights_[1]= weights_[13] = 0.0703660474881081247092674164507;
		weights_[2]= weights_[12] = 0.107159220467171935011869546686;
		weights_[3]= weights_[11] = 0.139570677926154314447804794511;
		weights_[4]= weights_[10] = 0.166269205816993933553200860481;
		weights_[5]= weights_[9] = 0.186161000015562211026800561866;
		weights_[6]= weights_[8] = 0.198431485327111576456118326444;
		weights_[7]= 0.202578241925561272880620199968;
		break;
	case 16:
		nodes_[0]=-0.989400934991649932596154173450; nodes_[15] = -1.0*nodes_[0];
		nodes_[1]=-0.944575023073232576077988415535; nodes_[14] = -1.0*nodes_[1];
		nodes_[2]=-0.865631202387831743880467897712; nodes_[13] = -1.0*nodes_[2];
		nodes_[3]=-0.755404408355003033895101194847; nodes_[12] = -1.0*nodes_[3];
		nodes_[4]=-0.617876244402643748446671764049; nodes_[11] = -1.0*nodes_[4];
		nodes_[5]=-0.458016777657227386342419442984; nodes_[10] = -1.0*nodes_[5];
		nodes_[6]=-0.281603550779258913230460501460; nodes_[9] = -1.0*nodes_[6];
		nodes_[7]=-0.0950125098376374401853193354250; nodes_[8] = -1.0*nodes_[7];
		weights_[0]= weights_[15] = 0.0271524594117540948517805724560;
		weights_[1]= weights_[14] = 0.0622535239386478928628438369944;
		weights_[2]= weights_[13] = 0.0951585116824927848099251076022;
		weights_[3]= weights_[12] = 0.124628971255533872052476282192;
		weights_[4]= weights_[11] = 0.149595988816576732081501730547;
		weights_[5]= weights_[10] = 0.169156519395002538189312079030;
		weights_[6]= weights_[9] = 0.182603415044923588866763667969;
		weights_[7]= weights_[8] = 0.189450610455068496285396723208;
		break;
	case 17:
		nodes_[0]=-0.990575475314417335675434019941; nodes_[16] = -1.0*nodes_[0];
		nodes_[1]=-0.950675521768767761222716957896; nodes_[15] = -1.0*nodes_[1];
		nodes_[2]=-0.880239153726985902122955694488; nodes_[14] = -1.0*nodes_[2];
		nodes_[3]=-0.781514003896801406925230055520; nodes_[13] = -1.0*nodes_[3];
		nodes_[4]=-0.657671159216690765850302216643; nodes_[12] = -1.0*nodes_[4];
		nodes_[5]=-0.512690537086476967886246568630; nodes_[11] = -1.0*nodes_[5];
		nodes_[6]=-0.351231763453876315297185517095; nodes_[10] = -1.0*nodes_[6];
		nodes_[7]=-0.178484181495847855850677493654; nodes_[9] = -1.0*nodes_[7];
		nodes_[8]= 0.0;
		weights_[0]= weights_[16] = 0.0241483028685479319601100262876;
		weights_[1]= weights_[15] = 0.0554595293739872011294401653582;
		weights_[2]= weights_[14] = 0.0850361483171791808835353701911;
		weights_[3]= weights_[13] = 0.111883847193403971094788385626;
		weights_[4]= weights_[12] = 0.135136368468525473286319981702;
		weights_[5]= weights_[11] = 0.154045761076810288081431594802;
		weights_[6]= weights_[10] = 0.168004102156450044509970663788;
		weights_[7]= weights_[9] = 0.176562705366992646325270990113;
		weights_[8]= 0.179446470356206525458265644262;
		break;
	case 18:
		nodes_[0]=-0.991565168420930946730016004706; nodes_[17] = -1.0*nodes_[0];
		nodes_[1]=-0.955823949571397755181195892930; nodes_[16] = -1.0*nodes_[1];
		nodes_[2]=-0.892602466497555739206060591127; nodes_[15] = -1.0*nodes_[2];
		nodes_[3]=-0.803704958972523115682417455015; nodes_[14] = -1.0*nodes_[3];
		nodes_[4]=-0.691687043060353207874891081289; nodes_[13] = -1.0*nodes_[4];
		nodes_[5]=-0.559770831073947534607871548525; nodes_[12] = -1.0*nodes_[5];
		nodes_[6]=-0.411751161462842646035931793833; nodes_[11] = -1.0*nodes_[6];
		nodes_[7]=-0.251886225691505509588972854878; nodes_[10] = -1.0*nodes_[7];
		nodes_[8]=-0.0847750130417353012422618529358; nodes_[9] = -1.0*nodes_[8];
		weights_[0]= weights_[17] = 0.0216160135264833103133427102665;
		weights_[1]= weights_[16] = 0.0497145488949697964533349462026;
		weights_[2]= weights_[15] = 0.0764257302548890565291296776166;
		weights_[3]= weights_[14] = 0.100942044106287165562813984925;
		weights_[4]= weights_[13] = 0.122555206711478460184519126800;
		weights_[5]= weights_[12] = 0.140642914670650651204731303752;
		weights_[6]= weights_[11] = 0.154684675126265244925418003836;
		weights_[7]= weights_[10] = 0.164276483745832722986053776466;
		weights_[8]= weights_[9] = 0.169142382963143591840656470135;
		break;
	case 19:
		nodes_[0]=-0.992406843843584403189017670253; nodes_[18] = -1.0*nodes_[0];
		nodes_[1]=-0.960208152134830030852778840688; nodes_[17] = -1.0*nodes_[1];
		nodes_[2]=-0.903155903614817901642660928532; nodes_[16] = -1.0*nodes_[2];
		nodes_[3]=-0.822714656537142824978922486713; nodes_[15] = -1.0*nodes_[3];
		nodes_[4]=-0.720966177335229378617095860824; nodes_[14] = -1.0*nodes_[4];
		nodes_[5]=-0.600545304661681023469638164946; nodes_[13] = -1.0*nodes_[5];
		nodes_[6]=-0.464570741375960945717267148104; nodes_[12] = -1.0*nodes_[6];
		nodes_[7]=-0.316564099963629831990117328850; nodes_[11] = -1.0*nodes_[7];
		nodes_[8]=-0.160358645640225375868096115741; nodes_[10] = -1.0*nodes_[8];
		nodes_[9]= 0.0;
		weights_[0]= weights_[18] = 0.0194617882297264770363120414644;
		weights_[1]= weights_[17] = 0.0448142267656996003328381574020;
		weights_[2]= weights_[16] = 0.0690445427376412265807082580060;
		weights_[3]= weights_[15] = 0.0914900216224499994644620941238;
		weights_[4]= weights_[14] = 0.111566645547333994716023901682;
		weights_[5]= weights_[13] = 0.128753962539336227675515784857;
		weights_[6]= weights_[12] = 0.142606702173606611775746109442;
		weights_[7]= weights_[11] = 0.152766042065859666778855400898;
		weights_[8]= weights_[10] = 0.158968843393954347649956439465;
		weights_[9]= 0.161054449848783695979163625321;
		break;
	case 20:
		nodes_[0]=-0.993128599185094924786122388471; nodes_[19] = -1.0*nodes_[0];
		nodes_[1]=-0.963971927277913791267666131197; nodes_[18] = -1.0*nodes_[1];
		nodes_[2]=-0.912234428251325905867752441203; nodes_[17] = -1.0*nodes_[2];
		nodes_[3]=-0.839116971822218823394529061702; nodes_[16] = -1.0*nodes_[3];
		nodes_[4]=-0.746331906460150792614305070356; nodes_[15] = -1.0*nodes_[4];
		nodes_[5]=-0.636053680726515025452836696226; nodes_[14] = -1.0*nodes_[5];
		nodes_[6]=-0.510867001950827098004364050955; nodes_[13] = -1.0*nodes_[6];
		nodes_[7]=-0.373706088715419560672548177025; nodes_[12] = -1.0*nodes_[7];
		nodes_[8]=-0.227785851141645078080496195369; nodes_[11] = -1.0*nodes_[8];
		nodes_[9]=-0.0765265211334973337546404093988; nodes_[10] = -1.0*nodes_[9];
		weights_[0]= weights_[19] = 0.0176140071391521183118619623519;
		weights_[1]= weights_[18] = 0.0406014298003869413310399522749;
		weights_[2]= weights_[17] = 0.0626720483341090635695065351870;
		weights_[3]= weights_[16] = 0.0832767415767047487247581432220;
		weights_[4]= weights_[15] = 0.101930119817240435036750135480;
		weights_[5]= weights_[14] = 0.118194531961518417312377377711;
		weights_[6]= weights_[13] = 0.131688638449176626898494499748;
		weights_[7]= weights_[12] = 0.142096109318382051329298325067;
		weights_[8]= weights_[11] = 0.149172986472603746787828737002;
		weights_[9]= weights_[10] = 0.152753387130725850698084331955;
		break;
	}
	return;
}

glni::quadrilateral::quadrilateral()
{
	initialized_ = false;
	size_ = -1;
}

glni::quadrilateral::quadrilateral(unsigned int size) : quadrilateral()
{
	initiate(size);
}

glni::quadrilateral::~quadrilateral()
{
	clear();
}

void glni::quadrilateral::clear()
{
	initialized_ = false;
	size_ = -1;
	g1d_.clear();
	return;
}

double glni::quadrilateral::integral(double xlow, double xhig, double ylow, double yhig, func_quad_ptr fx, void *att, void *instance)
{
#ifndef GLNI_NO_CHECK
	if (!initialized_)
		throw std::invalid_argument("glni::quadrilateral is not initialized.");

	if (xhig <= xlow || yhig <= ylow)
		throw std::invalid_argument("glni::quadrilateral::integral(...) got an invalid integration range.");
#endif

	double scl = (xhig-xlow)*(yhig-ylow)/4.0;
	double gx, gy, inte = 0.0;
	for (int i = 0; i < size_; i++)
	{
		for (int j = 0; j < size_; j++)
		{
			gx = 0.5*(xlow+xhig) + 0.5*(xhig-xlow)*g1d_.node(j);
			gy = 0.5*(ylow+yhig) + 0.5*(yhig-ylow)*g1d_.node(i);
			inte += g1d_.weight(i)*g1d_.weight(j)*scl*fx(gx, gy, xlow, xhig, ylow, yhig, att, instance);
		}
	}
	return inte;
}

void glni::quadrilateral::initiate(unsigned int size)
{
#ifndef GLNI_NO_CHECK
	if (size == 0 || size > 20)
		throw std::invalid_argument("glni::quadrilateral::initiate(...) got an invalid initiating order.");
#endif

	initialized_ = true;
	size_ = size;
	g1d_.initiate(size_);
	return;
}

glni::triangle::triangle()
{
	initialized_ = false;
	size_ = -1;
}

glni::triangle::triangle(unsigned int size) : triangle()
{
	initiate(size);
}

glni::triangle::~triangle()
{
	clear();
}

void glni::triangle::clear()
{
	initialized_ = false;
	size_ = -1;
	g1d_.clear();
	return;
}

double glni::triangle::integral(double x1, double x2, double x3, double y1, double y2, double y3, 
	func_tri_ptr fx, void *att, void *instance)
{
#ifndef GLNI_NO_CHECK
	if (!initialized_)
		throw std::invalid_argument("glni::triangle is not initialized.");

	if ((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1) <= GLNI_ZERO) // 这里没有设置为0是为了检测三个点是否共线
	{
		throw std::invalid_argument("glni::triangle::integral(...) got an invalid integration range.");
	}
#endif

	double ksi, eta;
	double gx, gy, inte = 0.0;
	double scl = x1*y2+x2*y3+x3*y1-x1*y3-x2*y1-x3*y2;
	for (int i = 0; i < size_; i++)
	{
		for (int j = 0; j < size_; j++)
		{
			ksi = 0.5*(1.0+g1d_.node(j));
			eta = 0.25*(1.0-g1d_.node(j))*(1.0+g1d_.node(i));
			gx = (1.0-ksi-eta)*x1 + ksi*x2 + eta*x3;
			gy = (1.0-ksi-eta)*y1 + ksi*y2 + eta*y3;
			inte += g1d_.weight(i)*g1d_.weight(j)*scl*fx(gx, gy, x1, x2, x3, y1, y2, y3, att, instance)
				*(1.0-g1d_.node(j))*0.125;
		}
	}
	return inte;
}

std::complex<double> glni::triangle::integral(double x1, double x2, double x3, double y1, double y2, double y3, 
	func_tri_ptr2 fx, void *att, void *instance)
{
#ifndef GLNI_NO_CHECK
	if (!initialized_)
		throw std::invalid_argument("glni::triangle is not initialized.");

	if ((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1) <= GLNI_ZERO) // 这里没有设置为0是为了检测三个点是否共线
	{
		throw std::invalid_argument("glni::triangle::integral(...) got an invalid integration range.");
	}
#endif

	double ksi, eta;
	double gx, gy;
	std::complex<double> inte = std::complex<double>(0.0, 0.0);
	double scl = x1*y2+x2*y3+x3*y1-x1*y3-x2*y1-x3*y2;
	for (int i = 0; i < size_; i++)
	{
		for (int j = 0; j < size_; j++)
		{
			ksi = 0.5*(1.0+g1d_.node(j));
			eta = 0.25*(1.0-g1d_.node(j))*(1.0+g1d_.node(i));
			gx = (1.0-ksi-eta)*x1 + ksi*x2 + eta*x3;
			gy = (1.0-ksi-eta)*y1 + ksi*y2 + eta*y3;
			inte += g1d_.weight(i)*g1d_.weight(j)*scl*fx(gx, gy, x1, x2, x3, y1, y2, y3, att, instance)
				*(1.0-g1d_.node(j))*0.125;
		}
	}
	return inte;
}

void glni::triangle::initiate(unsigned int size)
{
#ifndef GLNI_NO_CHECK
	if (size == 0 || size > 20)
		throw std::invalid_argument("glni::triangle::initiate(...) got an invalid initiating order.");
#endif

	initialized_ = true;
	size_ = size;
	g1d_.initiate(size_);
	return;
}

glni::cube::cube()
{
	initialized_ = false;
	size_ = -1;
}

glni::cube::cube(unsigned int size) : cube()
{
	initiate(size);
}

glni::cube::~cube()
{
	clear();
}

void glni::cube::clear()
{
	initialized_ = false;
	size_ = -1;
	g1d_.clear();
	return;
}

double glni::cube::integral(double x1, double x2, double y1, double y2, 
	double z1, double z2, func_cube_ptr fx, void *att, void *instance)
{
#ifndef GLNI_NO_CHECK
	if (!initialized_)
		throw std::invalid_argument("glni::cube is not initialized.");

	if (x2 <= x1 || y2 <= y1 || z2 <= z1)
		throw std::invalid_argument("glni::cube::integral(...) got an invalid integration range.");
#endif

	double gx, gy, gz, inte = 0.0;
	double scl = (x2-x1)*(y2-y1)*(z2-z1)/8.0;

	for (int i = 0; i < size_; i++)
	{
		for (int j = 0; j < size_; j++)
		{
			for (int k = 0; k < size_; k++)
			{
				gx = 0.5*(x1+x2) + 0.5*(x2-x1)*g1d_.node(k);
				gy = 0.5*(y1+y2) + 0.5*(y2-y1)*g1d_.node(j);
				gz = 0.5*(z1+z2) + 0.5*(z2-z1)*g1d_.node(i);
				inte += g1d_.weight(i)*g1d_.weight(j)*g1d_.weight(k)*scl
					*fx(gx, gy, gz, x1, x2, y1, y2, z1, z2, att, instance);
			}
		}
	}
	return inte;
}

void glni::cube::initiate(unsigned int size)
{
#ifndef GLNI_NO_CHECK
	if (size == 0 || size > 20)
		throw std::invalid_argument("glni::cube::initiate(...) got an invalid initiating order.");
#endif

	initialized_ = true;
	size_ = size;
	g1d_.initiate(size_);
	return;
}

glni::tetrahedron::tetrahedron()
{
	initialized_ = false;
	size_ = -1;
}

glni::tetrahedron::tetrahedron(unsigned int size) : tetrahedron()
{
	initiate(size);
}

glni::tetrahedron::~tetrahedron()
{
	clear();
}

void glni::tetrahedron::clear()
{
	initialized_ = false;
	size_ = -1;
	g1d_.clear();
	return;
}

double glni::tetrahedron::integral(double x1, double x2, double x3, double x4, 
	double y1, double y2, double y3, double y4, 
	double z1, double z2, double z3, double z4, 
	func_tet_ptr fx, void *att, void *instance)
{
#ifndef GLNI_NO_CHECK
	if (!initialized_)
		throw std::invalid_argument("glni::tetrahedron is not initialized.");
#endif // GLNI_NO_CHECK

	double ksi, eta, phi;
	double gx, gy, gz, inte = 0.0;
	double x12 = x2-x1, y12 = y2-y1, z12 = z2-z1;
	double x13 = x3-x1, y13 = y3-y1, z13 = z3-z1;
	double x14 = x4-x1, y14 = y4-y1, z14 = z4-z1;
	double crx = y12*z13-z12*y13;
	double cry = z12*x13-x12*z13;
	double crz = x12*y13-y12*x13;
	double scl = (x14*crx+y14*cry+z14*crz);
	if (scl < 0) scl = -1.0*scl; // 6 times volume

#ifndef GLNI_NO_CHECK
	if (scl < GLNI_ZERO)
		throw std::invalid_argument("glni::tetrahedron::integral(...) got an invalid integration range.");
#endif // GLNI_NO_CHECK

	for (int i = 0; i < size_; i++) // z
	{
		for (int j = 0; j < size_; j++) // y
		{
			for (int k = 0; k < size_; k++) // x
			{
				// evaluation point in the reference tetrahedron
				ksi = 0.125*(1.0+g1d_.node(k))*(1.0+g1d_.node(j))*(1.0+g1d_.node(i));
				eta = 0.125*(1.0+g1d_.node(k))*(1.0+g1d_.node(j))*(1.0-g1d_.node(i));
				phi = 0.25*(1.0+g1d_.node(k))*(1.0-g1d_.node(j));
				// map reference coordinates to the real ones
				gx = x12*ksi + x13*eta + x14*phi + x1;
				gy = y12*ksi + y13*eta + y14*phi + y1;
				gz = z12*ksi + z13*eta + z14*phi + z1;
				// Evaluate the integral in a reference tetrahedron, then multiple 6 times volume
				inte += g1d_.weight(i)*g1d_.weight(j)*g1d_.weight(k)*scl
					*fx(gx, gy, gz, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, att, instance)
					*(1.0+g1d_.node(k))*(1.0+g1d_.node(k))*(1.0+g1d_.node(j))/64.0;
			}
		}
	}
	return inte;
}

std::complex<double> glni::tetrahedron::integral(double x1, double x2, double x3, double x4, 
	double y1, double y2, double y3, double y4, 
	double z1, double z2, double z3, double z4, 
	func_tet_ptr2 fx, void *att, void *instance)
{
#ifndef GLNI_NO_CHECK
	if (!initialized_)
		throw std::invalid_argument("glni::tetrahedron is not initialized.");
#endif // GLNI_NO_CHECK

	double ksi, eta, phi;
	double gx, gy, gz;
	std::complex<double> inte = std::complex<double>(0.0, 0.0);
	double x12 = x2-x1, y12 = y2-y1, z12 = z2-z1;
	double x13 = x3-x1, y13 = y3-y1, z13 = z3-z1;
	double x14 = x4-x1, y14 = y4-y1, z14 = z4-z1;
	double crx = y12*z13-z12*y13;
	double cry = z12*x13-x12*z13;
	double crz = x12*y13-y12*x13;
	double scl = (x14*crx+y14*cry+z14*crz);
	if (scl < 0) scl = -1.0*scl;

#ifndef GLNI_NO_CHECK
	if (scl < GLNI_ZERO)
		throw std::invalid_argument("glni::tetrahedron::integral(...) got an invalid integration range.");
#endif // GLNI_NO_CHECK

	for (int i = 0; i < size_; i++)
	{
		for (int j = 0; j < size_; j++)
		{
			for (int k = 0; k < size_; k++)
			{
				ksi = 0.125*(1.0+g1d_.node(k))*(1.0+g1d_.node(j))*(1.0+g1d_.node(i));
				eta = 0.125*(1.0+g1d_.node(k))*(1.0+g1d_.node(j))*(1.0-g1d_.node(i));
				phi = 0.25*(1.0+g1d_.node(k))*(1.0-g1d_.node(j));
				gx = (x2-x1)*ksi + (x3-x1)*eta + (x4-x1)*phi + x1;
				gy = (y2-y1)*ksi + (y3-y1)*eta + (y4-y1)*phi + y1;
				gz = (z2-z1)*ksi + (z3-z1)*eta + (z4-z1)*phi + z1;
				inte += g1d_.weight(i)*g1d_.weight(j)*g1d_.weight(k)*scl
					*fx(gx, gy, gz, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, att, instance)
					*(1.0+g1d_.node(k))*(1.0+g1d_.node(k))*(1.0+g1d_.node(j))/64.0;
			}
		}
	}
	return inte;
}

void glni::tetrahedron::initiate(unsigned int size)
{
#ifndef GLNI_NO_CHECK
	if (size == 0 || size > 20)
		throw std::invalid_argument("glni::tetrahedron::initiate(...) got an invalid initiating order.");
#endif

	initialized_ = true;
	size_ = size;
	g1d_.initiate(size_);
	return;
}
